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In this paper we describe a method for evolving the projected Gross-Pitaevskii equation (PGPE) for a Bose 
gas in a harmonic oscillator potential. The central difficulty in solving this equation is the requirement that the 
classical field is restricted to a small set of prescribed modes that constitute the low energy classical region of 
the system. We present a scheme, using a Hermite-polynomial based spectral representation, that precisely im- 
plements this mode restriction and allows an efficient and accurate solution of the PGPE. We show equilibrium 
and non-equilibrium results from the application of the PGPE to an anisotropic trapped three-dimensional Bose 
gas. 
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I. INTRODUCTION 

Recently a variety of classical field methods have become 
popular in the description of ultra-cold Bose gases d El El El 
S|5| 011113 Mini 112 The 
appeal of these methods is that the dynamics of the modes 
are treated non-perturbatively so that non-equilibrium situa- 
tions (e.g. see 1 1 3 1) or strongly fluctuating equilibrium sys- 
tems (e.g. see 1 18 ]) can be accurately simulated. 

In Refs. 0Q2) we have developed a classical field theory, 
known as the Projected Gross Pitaevskii Equation (PGPE) for- 
malism, to describe the finite temperature Bose gas. This ap- 
proach has found good agreement with experiment in the crit- 
ical region of the condensation transition [18], and has seen 
numerous applications to regimes where traditional meanfield 
methods are inapplicable (e.g. see (jjSEQlEI]])- A key com- 
ponent of our theory (and the primary distinction from other 
finite temperature classical field theories |4|) that enables it 
to be applied to the quantitative description of experiments is 
the use of a projector, i.e. the explicit restriction of our de- 
scription to the low energy modes of the system. For typical 
regimes of interest of order a thousand modes of the system 
are sufficiently highly occupied to be treated using a classical 
field approach l22l . 

Over the past decade there has been extensive development 
of techniques for finding ground state solutions to the Gross- 
Pitaevskii equation and algorithms for evolving the conden- 
sate. The basic premise here is that the system is at zero tem- 
perature so that a single mode of the system is occupied and 
the underlying basis states used for the simulation are unim- 
portant as long as they span the condensate. As such, in gen- 
eral these methods are not immediately applicable to the finite 
temperature case. 

In Ref. lfT6l we have outlined a theoretical scheme for sim- 
ulating a finite temperature trapped Bose gas. In this paper 
we present our algorithm for that case as well as the uniform 
system, and discuss the techniques we use to initialize and an- 
alyze PGPE simulations. As the typical usage of the PGPE 
formalism is far removed from the well-known zero temper- 
ature Gross-Pitaevskii equation (used to describe pure con- 
densate dynamics), we use the results section to present a de- 
tailed example of application and analysis for the technique. 
This should be of benefit for others in the community wishing 



to adopt these techniques and also provides thermal quanti- 
ties for comparison that arise from well-characterized simu- 
lations. An important result we present is a demonstration of 
irreversible behavior of the trapped PGPE, giving evidence for 
re-thermalization of this system. 

The outline of this paper is as follows. In the remainder 
of this section we briefly introduce the PGPE evolution equa- 
tion. In Sec. [IT] we set up a convenient set of units and outline 
our generic spectral approach to solving the PGPE equation in 
a finite basis. The details of implementing the algorithm are 



first presented for a uniform system in Sec. Ill In this case 
the modes are plane-wave-like and the algorithm can be ef- 
ficiently implemented using Fourier transforms. In Sec. [IV 
the main result of the paper is presented: An implementation 
for the experimentally relevant case of a harmonically trapped 
system. Here the natural modes to work in are the harmonic 
oscillator eigenstates, and we show that a finite number of 
such modes can be propagated accurately and efficiently using 
appropriate quadrature grids to evaluate the nonlinear matrix 
elements. Results for the finite temperature evolution of a har- 
monically trapped system are presented in Sec. |V| before we 
conclude. 



A. PGPE theory 

The PGPE is a time-dependent nonlinear Schrodinger equa- 
tion of the form E E UM ' 



ih 



dt 



|wo|^| 2 ^}, 



where 



H sp = H + 6V(xL,t), 



H 



2m 



V 2 +F t rap(x). 



(1) 

(2) 
(3) 



H sp is the single particle Hamiltonian, which includes the 
dominant H part and we allow for a (small) perturbation part 
SV, ip = ^(x, t) is the classical matter wave field (taken to 
be normalized to unity), N is the total number of atoms de- 
scribed by the PGPE, Vtrap(x) is the external trapping poten- 
tial, and /To = ^irah 2 /m, with a the s-wave scattering length. 
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FIG. 1 : Schematic diagram showing the classical and incoherent re- 
gions of the single particle spectrum for a harmonically trapped Bose 
gas. The energy e cu t is usually chosen so that the average number of 
particles in the modes at the cutoff is n cu t 00 1. 



modes accurately. 

Most commonly used methods for propagating Schrodinger- 
type equations do not satisfy these requirements; in particular, 
many methods do not propagate all modes of the numerical 
basis faithfully. This leads to negligible errors if the highest 
modes are unoccupied, as is the case for the T = Gross- 
Pitaevskii equation. However, it is clear that methods based 
on such assumptions will not be appropriate for simulating 
fields where the dynamics of all the modes are important. 



B. Computational units 



The numerical implementation of the PGPE formalism poses 
a rather interesting challenge: Formally only the low energy 
modes of the system are classical (i.e. the classical region 
shown schematically in Fig. [T] l22l ) and should be the only 
modes retained in the numerical description, a restriction ex- 
pressed formally by the projector 

P{F(x)} = 2 0„(x) f d 3 x'C(x')F(x'), (4) 

where <j) n (x) are eigenstates of H and the summation is re- 
stricted to modes in the classical region. The action of V in 
Eq. ([4]) is thus to project the arbitrary function F(x) into the 
classical region (C), which we take to be defined by the sin- 
gle particle eigenstates up to some specified (single particle) 
energy e cut . For a full description of the system dynamics we 
will require a means of simulating the incoherent region (i.e. 
the modes complementary to the classical region - see Fig.[T]) 
and a means to couple the two regions. A theoretical formal- 
ism for doing this has been presented in Ref. fTOlL but has yet 
to be fully implemented numerically. However, the purpose of 
this paper is to present our approach for efficiently and accu- 
rately simulating Eq. ([I]) for the experimentally relevant case 
of a harmonic oscillator external potential. 

II. FORMAL ALGORITHM 

A. Numerical requirements 

The modes of the system are of central importance in the 
implementation of the PGPE and care must be taken in nu- 
merical implementations to ensure the modes are faithfully 
represented. Any useful simulation technique must satisfy the 
following requirements. 

(i) The space spanned by the modes of the simulation 
should match that of the classical region of the physi- 
cal system being simulated as closely as possible. That 
is, the modes should be the single-particle modes of the 
system (i.e. eigenstates of Hq) up to the prescribed en- 
ergy cutoff e cut . 

(ii) The assumption of high occupancy in all modes neces- 
sitates that the numerical scheme must propagate all 



For convenience we present the discussion of our numerical 
methods in computational units, indicated by tildes. We do 
this by introducing appropriate units of distance xq and time 
to in each of the following sections [41]. So for example, 
our dimensionless distance variable is defined as x = x/xo, 

dimensionless time is t = t/to, and classical field -0 = i/jx^ 2 . 
The coefficient of the nonlinear term in the Gross-Pitaevskii 
equation is given by the product NUq. In dimensionless units 
we define this as the nonlinearity constant C = NU^to/hx^. 

Adopting computational units, the projected Gross- 
Pitaevskii equation takes the form 

i 8 ^ = V{H SP ^ + C\^}. (5) 



C. Implementing the projector 

The projector can be implicitly implemented by restricting 
the classical field to the modes of interest, i.e. 

^(*, t) = ^ c n(t) </>n(x), (6) 

where {</> n (x)} are the eigenstates of #0 with respective 
eigenvalue e n . The projection is effected by limiting the sum- 
mation indices in Eq. ^ to the set of values 

C = {n : e n < e cut }, (7) 

i.e. the field ip only contains the modes of interest. 



1. Mode evolution 

Having restricted the modes for the purposes of the projec- 
tor, we can adopt these modes as our spectral basis and rep- 
resent their evolution exactly using a Galerkin approach (i.e. 
projecting Eq. ([5} on to our spectral basis). This leads to an 
evolution equation for the amplitudes 

dc 

= ~i [e n c n + F n + CG n ) , (8) 
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where 

F n = y"d 3 x^(x)jy(x,t)^(x,t), (9) 

G n = y"d 3 x0;(x)|^(x,t)| 2 ^(x,t), (10) 

are the matrix elements of the perturbation potential and non- 
linear term respectively. Once all the matrix elements on the 
right hand side of Eq. ^ are evaluated, the evolution of the 
system can be calculated using numerical algorithms for sys- 
tems of ordinary differential equations, e.g. the Runge-Kutta 
algorithm (e.g. see (23)). Here we do not concern ourselves 
with the particular choice of propagation algorithm, but in- 
stead focus on evaluating the matrix elements. We do mention 
in passing that by moving to an interaction picture, defined by 
the transformation c n (t) — exp(ie n t)c n (t), the explicit de- 
pendence on e n can be removed from Eq. ([8} (also see [ 24 1). 

III. IMPLEMENTATION FOR A UNIFORM SYSTEM 

Here we consider the numerical description of a Bose gas 
in a cuboid volume with linear dimensions {L x , L yj L z } and 
subject to periodic boundary conditions. To simplify our dis- 
cussion we consider the case where L x = L y = L z so that 
the ID basis states (see below) are identical in each direction. 
For a Bose gas in a uniform system with periodic boundary 
conditions the basis Hamiltonian takes the form 



with the boundary conditions, 



(2tt)2 



(11) 



$(x + l,y,z) = i/>(x,y + l,z) = t/j(x,y,z + l) = ij){x,y,z), 

(12) 

where we have taken xo = L x as the unit of length and to = 
mLl/irh as the unit of time. 



label the ID eigenstates, and note that the classically simu- 
lated region limits the values these can take to the set 



{a, /3, 7 : e a + ip + e 1 < e cut }, 



(17) 



i.e. a sphere of radius Ve cut in a/?7-space. For later conve- 
nience we define a max as the maximum value of a that oc- 
curs in C, i.e. the highest order basis state in each direction. 
For the planewave case we have a max ~ V^cut- This means 
that within the classical region there exists M = 2a max + 1 
distinct ID eigenstates (i.e. <p a ) in each direction (since 
a e [-cwx,cwx]), and M T « f M 3 3D basis states (0 n ) 
in the classical region. In what follows, we will adopt the 
notation ^ to indicate summations restricted to the classical 
region, i.e. all triplets of Greek indices in C. 



B. Evaluating the matrix elements 

The nonlinear matrix element ( [T0| takes the form 

G a(3l = J d 3 x^(x)^(y)^(z) ^(x,t) $(x,t). 

(18) 

Substituting the basis expansion for the field, given by Eqs. 
^ and (13) into Eq. ( 1_8) gives a series of matrix elements 
that can be independently integrated over the coordinate di- 
rections. Each of these ID integrals is of the form 



dx (p* a (x)<pg (x)(p v (x)<p a (x) , 



(19) 



which the indices in general can take any of the M values in 
the range -a max , -(a max - 1), ... , a max . These integrals 
are in some sense trivial for the planewave case, i.e. I a 5ua — 
5 a +s,v+a, however knowledge of these matrix elements is not 
sufficient to implement an efficient algorithm, as we discuss 
below. 



Expanding the field x/j in Eq. ([18| in terms of the basis 
states, the nonlinear matrix elements can be written as 



A. Separating into ID basis states 

The basis states (eigenstates of Ho) are separable into ID 
eigenstates, i.e. 



n (x) <-> <Pa(x)<pp(y)<p 7 (z), 



£n <"> Sa + 6p 
C n <— > Cqi/37 7 



where 



<Pa(x) 



(13) 
(14) 
(15) 



(16) 



are the ID eigenstates of the kinetic energy operator with the 
wavevectors k a chosen as harmonics of the perodicity inter- 
val, i.e. k a = 2ira, with a an integer, and have respective 
eigenvalues e a = a 2 . For clarity we use Greek subscripts to 



5Qr] v^p (Jtv 



(20) 



Explicitly carrying out the summations in Eq. ( 20 ) for all the 



matrix elements of G a p 1 requires 0(M 12 ) operations, and 
would be prohibitively slow for any practical calculation. We 
now show how a quadrature approach can be used to evaluate 
these integrals much more efficiently, requiring only 0(M 4 ) 
operations. The essence of this approach is to transform the 
field to a spatial representation where the nonlinear term is lo- 
cal. By choosing the spatial grid as an appropriate quadrature 
grid, the matrix elements will still be evaluated exactly. The 
basic idea of this approach (efficient evaluation of nonlinear 
matrix elements on a spatial grid) is widely used (e.g. see 
ES ESI [261 El EH El), although it is usually implemented 
in a manner that is not exact. This is typically an acceptable 
approximation if the highest energy modes of the system are 
unoccupied, a luxury not available in the PGPE. 
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In each spatial dimension, the quadrature grid of interest 
(for the uniform case) consists of Nq -points given by 



Xj=jAx, l<j<N Q 



(21) 



with spacing Ax — 1/Nq, which spans the spatial region 
(0,1]. The quadrature expression for an integral of an arbi- 
trary function / is 



J a 



Nq 

dxw(x) f(x) « ^2wjf(xj), (22) 



where w(x) is the weight function, and Wj are the quadrature 
weights. Formally such a quadrature can be implemented for 
our planewave case if we take a = 0, b = 1, w{x) = 1 and 
Wj = Ax. 

The requirement that our quadrature will exactly calculate 
the nonlinear matrix elements is equivalent to the requirement 
that the ID matrix elements fl9] ) are all evaluated exactly, i.e. 

Nq 

Iol^s = ^ Ax tp* (xj ) <p* {xj ) (p 1 (xj )<ps(xj), (23) 

= <W/3, 7 +5, (24) 

which holds for the quadrature described above if we take 
Nq > 2M. For what follows we take Nq = 2M. This 
quadrature can be applied to the 3D case, so that the nonlinear 
matrix elements are computed as 

ijk 



x\i/j(xi jk J)\ 2 ^(x ijk ,t), 



(25) 



where = yj,z k ), the y and z grids are taken to be 
identical to the x grid, and the summation over {i, j, k} here- 
after will be each taken to be over the range of values 1 to 
= 2M. 

Unless the perturbation potential can be expressed as a 
superposition of the planewave states ( [T6| ) of maximum 
wavevector /c max = 47ra max in each direction, the quadrature 
given above will not be exact for evaluating Eq. ([9]). However, 
since the perturbation will normally be small this approximate 
treatment should be satisfactory, i.e. we will take 

ijk 

x5V(5c ijk J)ijj(xi jk J). (26) 



C. Overview of numerical procedure 

Here we briefly overview how the quadrature described 
above can be efficiently implemented numerically. Starting 
from the basis set representation of the field (i.e. {c a/ 3 7 }) at 
an instant of time t, the steps for calculating the matrix ele- 
ments are as follows: 



1 . The field is transformed to a spatial representation ac- 
cording to 

tjj(5t ijk ,t) = ^ UjqUjpUkj Cg(3<y(t), (27) 

where the transformation matrices are defined as the ID 
basis states evaluated on the quadrature grid, i.e. 



Uia — (pa • 



(28) 



2. The integrands of the matrix elements ([9]) and ([10]) are 
then constructed, i.e. 

f(xijk) = SV(x ijk ,t)^(x ijk ,t), (29) 
g(iijk) = |^(x ij/c ,t)| 2 ^(x ijfe ,t). (30) 

3. Inverse transforming these integrand functions yields 
the desired matrix elements: 

F aPl = (Ax) 3 J2u* a U* U^f(5t ijk ), (31) 

ijk 

G aPl = (Ax) 3 J2uLU* U^g(i ijk ). (32) 

ijk 

The slowest step in this procedure is carrying out the basis 
transformation (steps 1 and 3). The computational cost of this 
step is 0(M 4 ) when carried out as a series of matrix multi- 
plications. However, for the planewave case this transforma- 
tion is equivalent to a fast Fourier transformation, which has 
a computational cost of O (M 3 log(M)) . In the last step the 
transforms to obtain F and G can be combined into a single 
transform, i.e. in practice we transform /(x^^) + Cg{S^ij k ), 
which yields F n + CG n (see Eq. ([5])). 



IV. IMPLEMENTATION FOR HARMONICALLY 
TRAPPED SYSTEM 

We now consider the main subject of this paper: The imple- 
mentation of the PGPE for the harmonically trapped system, 
i.e. where 

^trap(x) = X -m (uo 2 x x 2 + uo 2 y y 2 + lj 2 z z 2 ) . (33) 

In computational units the basis Hamiltonian, H , takes the 
form 



H 



1 M 1 

--V 2 + - 
2 2 



V 2 + ^(A^ 2 + A^ 2 + i 2 ), (34) 



where = co x /co z , X y = co y /co z , and we have used the 
harmonic oscillator frequency associated with the z-direction 
to define units of length xq = yjh/muj z and time to = u~ l . 
To simplify our discussion of the numerical method, we will 
take the harmonic trapping potential to be isotropic, i.e. \ x = 



1. This allows us to avoid using cumbersome notation 



to account for different spectral bases in each direction. 
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The same decomposition used for the planewave case ( [T3] ) 
can be applied in the harmonically trapped system if we 
take {<p a (x)} to be eigenstates of the ID harmonic oscillator 
Hamiltonian, i.e. 



~2dx 2 



1 



(35) 



with eigenvalue e a = (a + \), where we have taken a to 
be a non-negative integer. Such spectral representations have 
been considered previously for the zero-temperature (non- 
projected) Gross-Pitaevskii equation in Refs. l24l|29l . 

We define a max as the maximum value of a that occurs in 
C, i.e. for the harmonic case a max ~ e cut — \ . This means 
that within the classical region there exists M = a ma x + 1 
distinct ID eigenstates (i.e. (p a ) in each direction (since a 
takes the values 0, 1, ... , a max ), and M T « ^M 3 3D basis 
states (0 n ) in the classical region. 



A. Oscillator state properties 

We briefly review the properties of the harmonic oscillator 
states, which we utilize as our spectral basis. 



The eigenstates of the ID single particle Hamiltonian (35) 



are 



<p a (x) = h a H a (x)e~ 



2 /2 



(36) 



where h a — [2 a a!v / 7r] -1 / 2 is the normalization constant, and 
H a (x) is a Hermite polynomial of degree a, defined by the 
recurrence relation 

H a +i(x) = 2xH a (x) - 2aH a -i(x), a = 1,2,... (37) 

with Ho(x) = 1, and H\{x) = 2x. The harmonic oscillator 
states are eigenstates of the Fourier transform operator with 
eigenvalue (— i) a , i.e. 



<Pa(P) = (-*)" 



dx ( 



C (fa(x). 



(38) 



Thus knowledge of the basis amplitudes c a/ 3 7 allows us to effi- 
ciently and precisely construct the momentum representation 
of the classical field, i.e. 

*(P) = E(-0 a+/3+7 C«/3 7 Pa(p*)^(p„)£y(Pz). (39) 

a (3 j 



1. Step operators 

It is useful to consider the so-called step operators of quan- 
tum mechanics, defined as 



(- 



d_ 

V2 \ dx 
1 / d 



V2 \dx 



(40) 
(41) 



which are mutually adjoint and have the commutation relation 

[a-,a+] = 1. (42) 

In a similar manner we can define step operators for the y and 
z directions, although for the sake of brevity we avoid doing 
this here. 

Applying the step operators to the ID eigenstates yields 



a+(p a (x) = y/a + l(p a+1 (x), 
a~(p a (x) = y/a<p a -i(x), 



(43) 
(44) 



so that the matrix representation of these operators in the spec- 
tral basis is 

{at) a(3 = J dx<pl(x)a+<pi3(x), (45) 

= Vi9+T*a,/3+i, (46) 
{K) a p = V^^-i. (47) 

Most importantly this allows us to represent the operators x 
and d x ~ d 



g~ in the spectral basis (exactly) as 

(x)clP = A=(a>5 + a>x)*Pi 



V2 



13 + 1 



"^a,/3-lj 



(d x ) 

V J at 



1 



(3 + 1 



/3+1- 



(48) 

(49) 
(50) 

(51) 



So, for example, consider the position expectation of the field, 
i.e. (x(t)) = j <i 3 xx|^(x, t)| 2 . Note that this expectation is 
a quantum mechanical average at time t, rather than a thermal 
(ensemble/time) average. This quantity can be calculated in 
the spectral basis as 



(x(t)) 



^2 c hi 



(*) 



it), 



(52) 



where indicates a restricted summation over the vari- 

ables {a, $, 7, 5} such that both {a, /3, 7} and {(5, (3, 7} lie in 
C. While this appears to be of computational cost 0(M 4 ), in 
fact the sparseness of the underlying step operators means that 
this operation is 0(M 3 ). 

Some care needs to be taken when applying the step op- 
erators. Formally, these operators have infinite size matrix 
representations, whereas our classical field application is in- 
trinsically finite by virtue of the energy cutoff. An implication 
of this is that, e.g. the action of the a+ operator will take the 
highest x-modes in the classical region onto modes above the 
cutoff, which would then be lost from our description. This 
is a problem when we wish to consider the product of oper- 
ators and can be avoided by using the commutation relation 



( 42 ) to write operators of interest in a normally ordered form, 
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whereby lowering operators (a x ) occur before raising opera- 
tors (a+). For example, 

^ = \{at + a~)\ (53) 
= ll(at) 2 + (a-) 2 + 2aia- + l], (54) 

where we have replaced a~a+ by a+a~ + 1 on the second 
line according to our normal ordering prescription. 

We do not develop these ideas any further here, but empha- 
size the utility of these results in the analysis of classical field 
calculations. Many common operators, and hence matrix el- 
ements or expectations of interest, can be expressed in terms 
of products of position and momentum (derivative) operators 

(e.g. angular momentum L z = —i ^xd y — yd x ^). These can 

be reformulated in terms of the normally ordered step opera- 
tors that can be applied to the spectral basis as 0(M 3 ) oper- 
ations, and are exact. This avoids the 0(M A ) cost of trans- 
forming to quadrature grids and the associated difficulties of 
coming up with accurate expressions on such grids for the op- 
erators of interest. 



usual weight function for Gauss-Hermite quadrature, the in- 
tegral can be exactly evaluated using a three-dimensional spa- 
tial grid of 8 (M - l) 3 points (i.e. 2 (M — 1) points in each 
direction (42J), i.e. 

{xi,Xj,Xfc), (59) 

ijk 

where X{ and Wi are the 2 (M — 1) roots and weights of the 
ID Gauss-Hermite quadrature with weight function w(x) = 
exp(-2x 2 ) (30). 

The perturbation potential can be calculated exactly on this 
grid if it is of the form 

5V(5t, t) = e~^ 2+ y 2+s2) R(x, y, z), (60) 

where R(x, y, z) is a polynomial of maximum degree 2(M — 
1) in the independent variables. However, if we assume the 
perturbation is small, then it will be permissible to evaluate it 
approximately on the same quadrature grid used for the non- 
linear term. 



C. Overview of numerical procedure 



B. Evaluating the matrix elements 



The nonlinear matrix element (10) takes the form given in 
Eq. ( [T8] ). An important observation made in Ref. l29l is that 
these matrix elements can be computed exactly with an appro- 
priately chosen Gauss-Hermite quadrature. To show this we 
note that because the harmonic oscillator states are of the form 
— h a H a (x) exp(— x 2 /2), where H a (x) is a Hermite 
polynomial of degree a, the field (at any instant of time) can 
be written as 



^(x,t) = Q(x,y,z)e~ 



(x 2 +y 2 +z 2 )/2 



(55) 



where 



Q(x,y,z) = ^2c a(3l (t) h a H a (x)h{3H{3(y)hiH^(z), 

(56) 

is a polynomial that, as a result of the cutoff, is of maximum 
degree M — 1 in the independent variables. 

Similarly, it follows that because the interaction term flO} 
is fourth order in the field, it can be written in the form 

G a(3l = J d 3 xe- 2 ^+y 2 +~ z ^P a ^(x,y,z), (57) 

where 

P a(3l (x,y,z) = h a H a (x)hpHi3(y)h 7 H 7 (z) 

x\Q(x,y,z)\ 2 Q(x,y,z), (58) 

is a polynomial of maximum degree 4 (M — 1) in the inde- 
pendent variables. Identifying the exponential term as the 



Here we briefly overview how the quadrature described 
above can be efficiently implemented numerically. This dif- 
fers slightly from the planewave case by virtue of the less 
trivial nature of the quadrature weights and weight functions. 
Starting from the basis set representation of the field (i.e. 
{cq;/3 7 }) at an instant of time t, the steps for calculating the 
matrix elements are as follows: 

1. We transform the field to a spatial representation ac- 
cording to 

xjj(5t ijk ,t) = ^ U i a UjpU k >y c a(3^(t), (61) 
a/37 

where the transformation matrices are defined as the ID 
basis states evaluated on the quadrature grid, i.e. 



(62) 



2. The integrands of the matrix elements ([9]) and ( fT0| ) are 
then constructed as quadratures by appropriately divid- 
ing by the weight function and pre-multiplying by the 
weights l43lL i.e. 

f(xi jk ) = WiWjW k e 2 ^ kl ' ''5V(x ijk ,t)ijj(x ijk ,t), (63) 
g(xi jk ) = WiWjW k e 2 ^ k \ 2 \^ ijk ,i)\ 2 

3. Inverse transforming these integrand functions yields 
the desired matrix elements: 

Fa^ = Y. U ioPh U k 1 f^ijk), (65) 

ijk 

Ga^ = Y. U ioPh U k 1 9{^ijk). (66) 

ijk 
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The slowest step in this procedure is carrying out the basis 
transformation. The computational cost of this step is 0(M 4 ) 
when carried out as a series of matrix multiplications. 



As shown in Eq. ( [39] ) in the spectral basis there is a sim- 
ple relationship between the position and momentum repre- 
sentations. Thus to obtain the momentum space field we use 
Ufa — (—i) a( Pa(Pi) as the transformation matrices, i.e. 



D. Other transforms 

For completeness we describe how the basis transforma- 
tions can be generalized for use in preparing states in the spec- 
tral basis and transforming to position and momentum grids. 
While being unrelated to propagation, these procedures are 
part of an essential set of tools for the general application of 
the method and analysis of results. 



1. Projecting a position space state onto the spectral basis 

Frequently we are presented with a field specified in the 
position basis, i.e. ^ A (x), and need to obtain its spectral rep- 
resentation, i.e. 



(67) 



If the modes of the classical region span the function ip A (x) 
then this procedure provides an exact representation, however, 
in general the resulting coefficients represent the func- 
tion projected into the classical region, i.e., V{^ A }. 

In practice we implement this transform onto the spectral 
basis in a similar manner to that used in Eq. ([65]), i.e., as 



ijk 



(68) 



where the transformation matrices are defined as the ID basis 
states evaluated on the quadrature grid, i.e., Ui a = <p a (%i) 
as in Eq. ([62] but now the quadrature is different to that used 
to evaluate the nonlinear matrix elements. Replacing ijj A by 
V{^ A } in Eq. (67 ) we note that integrand is of the form of a 
exponential exp(^|x| 2 ) times a polynomial function of max- 
imum degree 2 (M — 1) in each of the coordinates. These 
integrals are exactly computed in each direction using Gauss- 
Hermite quadrature defined for the weight function w(x) = 
exp(— x 2 ) with M roots Xi and weights Wi. 



2. Transforming the classical field to arbitrary position and 
momentum grids 

An important aspect of any classical field method is the 
ability to transform the results to desired grids for analysis 
and/or visualization. This is particularly convenient for the 
spectral representation we have adopted here. The basic pro- 
cedure is as indicated in Eq. ( [6T) except that the grid points are 
now arbitrary and need not be related to any quadrature grid. 
Taking one direction to have a single point, e.g., {xj} = x\, 
allows us to take slices of the classical field. 



U! a ULU p k c a0 ,(t). 



(69) 



3. Column densities 

Computing column densities is important for modeling 
ultra-cold Bose gases, since this is what is measured in ex- 
periments using absorption imaging (the primary method of 
analyzing these systems). For example, when the in situ sys- 
tem is imaged using light propagating along the z-direction, 
the observable corresponds to the column density defined as 



dz|^(x)| 2 



(70) 



On the other hand, often the trapping potential is turned off 
and the system is allowed to expand freely before it is im- 
aged. If the expansion time is sufficiently long the measured 
signal is related to a momentum space column density, e.g. 
for imaging along z 



n c (PxiPy) 



I 



dp z |$(p 



~\|2 



(71) 



If the imaging direction (i.e. direction along which we inte- 
grate) corresponds to a coordinate direction the column den- 
sity can be conveniently and exactly evaluated. We will con- 
sider the case in Eq. ( [70} for definiteness, though the same 
procedure immediately applies to the momentum case, and 
for column densities taken along other axes. 

We transform the field to a position grid ^(x^fe), where the 
x and y grids can be chosen arbitrarily, but the z grid needs 
to be a M point quadrature grid of the type discussed below 



Eq. (68 ). This choice for the z grid ensures that 



h c (xi,yj) = ^ w fc e^|?/i(x ijfc )| s 



(72) 



exactly computes the z integral in Eq. ( 70 ) where the weights 



are those described below Eq. ( 68 ). 



E. Three body term 

Another important term that often needs to be evaluated 
(e.g. see Ref. |31|) is the so called three-body term, which 
can be implemented by adding the term — ^iKs^\ 4 ip to the 
right hand side of Eq. where K 3 characterizes the loss of 
atoms due to three-body recombination. Generally this pro- 
cess is small in comparison to the two-body interaction term 
(of strength C), and causes the field to lose normalization. 
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Including this term in the evolution equation ([5]) we obtain 

dc 1 ~ 

= ~i [e n c n + F n + CG n ] - -K 3 Jn, (73) 

where we have introduced the matrix element 

Jn = J d 3 xC(x)|^(x,t)|^(x,f). (74) 

Decomposing this into ID eigenstate basis (i.e. n — > {a/37}) 
we have that 

Jc/3 7 = J d 3 xe- 3 ^+y 2 +^S af3l (x,y,z), (75) 

where S a p>y is a polynomial of maximum degree 6(M — 1) in 
the independent variables. 

Using the same arguments as applied to the normal inter- 
action term (see Sec. |IIIB| ), we can show that this term could 
be evaluated exactly on a three-dimensional spatial grid con- 
sisting of 27 (M — l) 3 points (i.e 3(M — 1) points along each 
direction) using a quadrature appropriate to the weight func- 
tion w(x ) = exp(— 3x 2 ) in each direction. 

However, like the perturbative potential, usually this term 
has a rather small effect on the dynamics and the additional 
expense of implementing a special transform to exactly eval- 
uate the J-matrix elements is unnecessary. In such cases it 
should be an acceptable approximation to calculate these ma- 
trix elements on the same grids used for the two-body interac- 
tion. 



F. Convergence and future numerical work 



Relative Tolerance Number of steps 


SN 




SX 




1 x 10" 2 


131 


2.54 x 10" 


-2 


3.13 x 10" 


-1 


5 x 10" 3 


135 


2.04 x 10" 


-2 


1.94 x 10" 


-1 


1 x 10" 3 


186 


6.72 x 10" 


-3 


2.01 x 10" 


-2 


5 x 1(T 4 


213 


3.52 x 10" 


-3 


5.39 x 10" 


-3 


1 x 10" 4 


294 


7.13 x 10" 


-4 


2.25 x 10" 


-4 


5 x 10" 5 


335 


3.57 x 10" 


-4 


5.54 x 10" 


-5 


1 x 10" 5 


455 


7.02 x 10" 


-5 


2.15 x 10" 


-6 


5 x 10" 6 


527 


3.50 x 10" 


-5 


5.29 x 10" 


-7 


1 x 10" 6 


725 


6.93 x 10" 


-6 


2.05 x 10" 


-8 


5 x 1(T 7 


829 


3.46 x 10" 


-6 


4.98 x 10" 


-9 


1 x 10" 7 


1139 


6.92 x 10" 


-7 


1.67 x 10" 


-10 


5 x 10" 8 


1311 


3.45 x 10" 


-7 


3.28 x 10" 


-11 


1 x 10" 8 


1797 


6.92 x 10" 


-8 


N/A 





TABLE I: Convergence properties of evolution algorithm. The rela- 
tive error tolerance, number of steps needed to obtain that error toler- 
ance, and the final error measures SN and SX are shown (see text). 
Simulation with 10 ~ 8 tolerance is used as cf in measuring error of 
less accurate calculations. Other parameters: t = 2n, C = 1000, 
X x = 4, X y — 1, e C ut = 35 and the initial state has energy E = 18 
(see Sec.0. 



Here we present some evolution convergence results for 
our algorithm and discuss avenues for future development 
of PGPE numerics. Unlike the zero temperature Gross- 
Pitaevskii equation, where convergence can be ensured by in- 
creasing the basis size, here such a procedure will increase 
the classical region size (number of degrees of freedom) and 
hence effect the thermal properties of the system (44). Thus 
for the PGPE formalism the essential requirement is that all 
modes within the classical region are propagated precisely. 

We have used an adaptive step Runge-Kutta-Fehlberg algo- 
rithm to evolve the classical field equation ([3]) with a speci- 
fied relative error tolerance. Since computing the matrix ele- 
ments for our harmonically trapped algorithm is of computa- 
tional cost 0(M 4 ) the development of higher order or more 
efficient propagation algorithms would be desirable (e.g. see 
Refs. (32j [33] |34j |35) for various algorithms developed for 
evolving the zero temperature Gross-Pitaevskii equation). For 
example, the simulation presented in column (d) of Figs. [2]|4] 
has a basis of Mt = 1857 modes propagated over a time in- 
terval of t ~ 2007T. At a relative tolerance of 10 ~ 6 this simu- 
lation took approximately 2hrs and 23 minutes on a 2.66GHz 
Xeon (Woodcrest) processor, consisting of 61,147 steps in- 
cluding 3,013 failed steps. In Table [i] we examine the evolu- 
tion algorithm convergence properties using the two measures 

M T 

SX = Y,\ Cj (t)-cf(t)\ 2 : (76) 

M T 

SN = l-g| Cj -(f)| 2 , (77) 

i.e. a difference measure of the final states and the loss in nor- 
malization, where cf(t) are the mode amplitudes at time t of 
a more accurate simulation. 

While we have strived to ensure highly accurate time propa- 
gation, statistical mechanical arguments suggest that this may 
be unnecessary and some additional noise from imprecise evo- 
lution may assist the system explore the ensemble of available 
states more rapidly, as long as this noise does not effect the 
constants of motion describing the system. For our evolution 
algorithm the normalization of the field is not conserved and 
anecdotal evidence suggests that monitoring changes in nor- 
malization can be used as a summative assessment of the evo- 
lution accuracy, e.g. for the simulation described above the 
final state normalization was 0.99989. 

Here we have considered a spectral cutoff in the ideal sin- 
gle particle basis, i.e. harmonic oscillator states. (We refer 
readers to Ref. 1 36 ] for a discussion of how numerics can be 
implemented for a rotating trapped system). For sufficiently 
high energy cutoff, this basis approximately diagonalizes the 
many-body problem. However, for situations requiring lower 
energy cutoffs (e.g. when T <C T c ), the Hartree-Fock inter- 
acting modes would form a better basis for defining the clas- 
sical region. Since these modes do not have an associated 
quadrature, the scheme presented here cannot be immediately 
extended to this regime. Another pathway, which may allow 
more flexibility in the nature of the cutoff, would be to use 
a more efficient representation (e.g. plane-wave/Fourier trans- 
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form) and approximately implement the projection. 

V. RESULTS 
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FIG. 2: (color online) Position and momentum density slices in the 
x2/-plane. Row (al)-(dl) instantaneous position density slices for 
z = 0. Row (a2)-(d2) time-averaged position density slices for z — 
0. Row (a3)-(d3) instantaneous momentum density slices for p z = 0. 
Row (a4)-(d4) time-averaged momentum density slices for p z = 0. 
Parameters: Column (al)-(a4) E = 11.64, / con d = 0.87, column 
(bl)-(M) E = 15.07, /cond = 0.61, column (cl)-(c4) E = 20.70, 
/ cond = 0.24, column (dl)-(d4) E = 25.83, / con d = 0.01. Other 
parameters: C — 750, X x = 4, X y = 1 and e cu t — 35. 

In this section we consider the application of the PGPE for- 
malism to a harmonically trapped system with a fixed energy 
cutoff, and take SV = 0. To elucidate features of our method 
we consider an anisotropic harmonic trap with X x — 4 and 
X y = 1, i.e. the system has a fat pancake geometry with the 
x-direction being tightly confined. Taking 6 cu t = 35 there 
are M T = 1857 single particle modes in the classical region 
with the maximum order of eigenstate occurring in each co- 
ordinate direction being (a max )^ = 8, (a max ) y = 32, and 

(^max)^ = 32. 

The PGPE is ergodic and the generic method to study fi- 
nite temperature regimes is to begin with a randomized initial 
state with some definite energy, E, as quantified by the energy 
functional 

E$] = J d 3 x^*|^ sp + ^|^| 2 J^. (78) 

This energy is a constant of motion for the PGPE ^ and 
forms a convenient macroscopic constraint for specifying the 
thermal state of the system. The procedure for making such 
energy states is rather arbitrary. We choose to make use of 
the Thomas-Fermi approximation to the ground state of the 
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FIG. 3: (color online) Position and momentum density slices in the 
2/z-plane. Row (al)-(dl) instantaneous position density slices for 
x = 0. Row (a2)-(d2) time-averaged position density slices for x = 
0. Row (a3)-(d3) instantaneous momentum density slices for p x = 0. 
Row (a4)-(d4) time-averaged momentum density slices for p x = 0. 
Parameters: Column (al)-(a4) E = 11.64, / con d = 0.87, column 
(bl)-(b4) E = 15.07, / con d = 0.61, column (cl)-(c4) E = 20.70, 
jf cond = 0.24, column (dl)-(d4) E = 25.83, / con d = 0.01. Other 
parameters: C — 750, X x = 4, X y = 1 and e cu t — 35. 
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FIG. 4: (color online) Position and momentum column densities in- 
tegrated along the 5-axis. Row (al)-(dl) instantaneous position col- 
umn densities. Row (a2)-(d2) time-averaged position column den- 
sities. Row (a3)-(d3) instantaneous momentum column densities. 
Row (a4)-(d4) time-averaged momentum column densities. Param- 
eters: Column (al)-(a4) E = 11.64, / con d = 0.87, column (bl)- 
(M) E = 15.07, /cond = 0.61, column (cl)-(c4) E = 20.70, 
/ cond = 0.24, column (dl)-(d4) E = 25.83, / con d = 0.01. Other 
parameters: C — 750, X x = 4, X y = 1 and e cu t — 35. 
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energy functional, given by 



^TF(X) = y A*TF-Wx) g ( MTp _ y tiap{ ^ (?9) 

where is the unit step function and //tf = 

^_(15A :E A 2y C/47r) 2 / 5 is the Thomas-Fermi chemical potential 
1371. We mix the Thomas-Fermi state (projected onto the 



spectral basis according to the procedure given in Sec. IV D 1 ) 
with a high energy randomized state B31 in the appropriate 
ratio to obtain a normalized state of the desired energy. Such 
states, for a range of energies, are propagated according to the 
PGPE for a time period of t = 200tt and 500 states of the field 
are saved at equally spaced times during the evolution, i.e. at 
times tj = 2tt7'/5, with < j < 500. 

Due to the stochastic nature of the initial condition for the 
simulation the instantaneous behavior of the system is of re- 
duced importance, instead the quantities of interest are the 
macroscopic observables that can be computed from the sys- 
tem evolution. We consider a few such observables here, and 
refer the reader to Ref. [16] for a more general discussion on 
this topic. 

One observable we will concern ourselves with is the aver- 
age position density calculated as the time-average 



Wx)| 2 



x,*;)l' 



(80) 



where the summation is taken over some subset of M s saved 
states. Here we use the last 350 of the saved states for averag- 
ing. We use the first part of the evolution (i.e. up to t = 607r) 
to allow the arbitrarily chosen initial state to relax to equilib- 
rium, a point we consider further in Sec. |VB| Similarly we 
can construct the momentum density by time- averaging the 
Fourier transformed classical field ( [39] ). 

Another quantity of interest is the condensate fraction. The 
definition we use was provided by Penrose and Onsager (38), 
and identifies the condensate fraction f COYl d as the largest 
eigenvalue of the one-body density matrix, defined in terms 
of the field as G ljB (x, x') = (?/i*(x)^(x / )). In our formalism 
this quantity is equivalently and much more efficiently com- 
puted in the mode basis as 



G 



IB 



( C m C n) 5 



(81) 



where () is evaluated using time-averaging [e.g., see Eq.[80). 
It is also possible to calculate thermal parameters from the dy- 
namical evolution of the field, such as temperature and chem- 
ical potential, and we refer the reader to Ref. ifTTl for more 
details. 



A. Instantaneous and time averaged density profiles 

Our results for the instantaneous and averaged density pro- 
files are shown in Figs.[2]and[3] In Fig.|2]the system is viewed 
in the xy-plane using a z = slice through the field. The posi- 
tion space density clearly reveals the anisotropy of the confin- 
ing potential. As the energy of the system increases it fluctu- 
ates more strongly (as revealed in the instantaneous density 



slices), however the averaged density profiles change quite 
gradually. In contrast, the momentum space density (shown 
as p z =0 slices) changes much more significantly with sys- 
tem energy. At low energies [e.g., Fig. |5Ja3) and (a4)] the 
condensate is prominent in the system, and is easily identi- 
fied as the dense anisotropic momentum peak centered at zero 
momentum. This momentum anisotropy is the opposite of 
that observed in position space due to the Heisenberg rela- 
tionship between position and momentum size of the conden- 
sate wavepacket along each direction. We also see a back- 
ground radially symmetric distribution in the averaged mo- 
mentum density [e.g., Fig.[2ja4)] that increases in prominence 
(at the expense of the condensate) as the energy of the system 
increases [e.g., see figures along the row Fig.|2ja4)-(d4)]. This 
background constitutes what is usually referred to as the ther- 
mal cloud of the system - or at least the component of this 
that exists within the classical region [46 ] . Comparing the av- 
eraged and respective instantaneous momentum density pro- 
files [i.e., see Fig. [2] along the rows (instantaneous) (a3)-(d3) 
and (averaged) (a4)-(d4)] we see that the modes contributing 
to this radially symmetric density fluctuate quite strongly, and 
adequate time-averaging to obtain their averaged profiles is 
essential. 

In Fig.[3]the system is viewed in the ?/z-plane using ax = 
slice through the field (p x =0 slice for the momentum space 
results). The averaged position and momentum space density 
are isotropic like the trapping potential in this plane. These 
results also allow us to appreciate the role of the classical 
region cutoff on typical fluctuations in the system. Com- 
pare the instantaneous momentum densities in Fig.|2jd3) and 
Fig. [3jd3). In Fig. |3jd3) the length scale for fluctuations is 
the same in the p y and p z directions, reflecting the symmetry 
of the potential. In Fig.|2jd3) the length scale for fluctuations 
is much longer in the p x direction than in the p y direction. 
We do not see the same features in the position space fluc- 
tuations, e.g. in Fig. [2jdl) and Fig. |3jdl) where the length 
scales in all directions appear to be similar. These observa- 
tions can be explained by a simple argument based on a semi- 



classical analysis of the single particle Hamiltonian (34): The 
maximum position extent of the system in the j -direction is 
Xj ~ \/2e cu t/Aj, and using that the number of modes in this 
direction is Mj ~ e cut / Xj , we obtain that the typical fluctu- 
ation length scale is Axj ~ ^/2/ e cut l47lL i.e., independent 
of trap frequency (i.e. Xj) in position space. In momentum 
space we instead have that the momentum extent of the field, 
Pj ~ V2e cu t, is independent of the trap frequency, while 
the length scale for fluctuations Apj ~ Xj ^2/ e cut is pro- 
portional to the trap frequency. This qualitatively agrees with 
the observation in Fig. |2jd3) that the momentum fluctuation 
length scale is about 4 times longer in the p x direction than in 
the p y direction. 

For comparison with the density slices we have just dis- 
cussed in Fig. [4] we show the instantaneous and time-averaged 
column densities in position and momentum space. Taking 
the column density in some sense averages the fluctuations in 
the instantaneous system state compared to the density slice. 
The instantaneous column densities closely approximate the 
images taken in experiment, and thus the fluctuations seen 
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should be measurable, however we note that our results here 
exclude the contribution from the incoherent region which in 
general will be important to include. 



B. Relaxation to equilibrium 

Finally we consider some dynamical aspects of how a dis- 
torted system, described by the PGPE, returns to an equilib- 
rium state. To do this we construct two initial states with 
the same total energy but rather different density distributions. 
Density slices of these two states are shown in Figs. |5ja) and 
(c). Our preparation procedure for those two states is as de- 
scribed before, except that here we choose to use a Thomas- 
Fermi state that is distorted so as to be extended along the 
x and y directions, respectively. As conjectured earlier, since 
both states have the same energy they should relax to the same 
(dynamical) equilibrium. In Figs. |5Jb) and (d) we show the 
density slices of these states after propagating the system for 
30 (z-direction) trap periods. The density slices of these final 
states qualitatively look much more similar to each other than 
the initial states did, consistent with both systems reaching the 
same equilibrium. However, we can provide other evidence 
for the system returning to equilibrium. In Figs.[5je) and (f) 
we examine the position variance of the classical field along 
the y direction defined as 



■(y(t)) = (y(t) 2 >-(y(t)> 2 



(82) 



i.e. the quantum mechanical average of the field at time t. This 
provides a measure of the system width in that direction. We 
note that this quantity is efficiently evaluated with the step op- 
erator formalism discussed in Sec. |IV A 1 Fig. [5je) reveals 
that var (y(t)) initially shows strong evolution that is charac- 
teristic of the initial state being far from equilibrium. This ini- 
tial evolution is strikingly different for the two cases reflecting 
that initial states considered here are respectively compressed 
[Fig. |5ja)] and extended [Fig. [5jc)] in the y direction. After 
a few trap periods the dynamics of var(y) is seen to reduce 
considerably. 

In Fig. [5jf) we examine dynamics of var(y) after the system 
has relaxed for 10 trap periods. The fluctuations in var(y) are 
much reduced in amplitude compared to what was observed 
in the initial evolution. Indeed, the fluctuations here are con- 
sistent with thermal fluctuations of the equilibrium state, and 
it is clear that both initial states considered have relaxed to an 
equilibrium state with approximately the same average value 
for var(y). We also note that both systems end up with ap- 
proximately the same condensate fraction of f con d ~ 0.15, 
obtained from density matrix calculated by time-averaging 
states over the last 10 trap periods of the system evolution. 



VI. CONCLUSIONS AND OUTLOOK 

In this paper we have described an efficient spectral method 
for solving the projected Gross-Pitaevskii Equation for a Bose 
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FIG. 5: (color online) Relaxation to equilibrium. Density slices of 
the two non-equilibrium initial states (a) and (c), and the respective 
states they evolve to at t = 60tt (b) and (d). Both states have E — 
22.00. The variance of the classical field in the ^-direction (e) during 
the first few trap periods and (f) after 10 trap periods. Result for 
initial state (a) (black line) and initial state (b) (grey line). Other 
parameters: C — 750, \ x — 4, X y = 1 and e cu t — 35. 



gas in a harmonic oscillator potential. We have shown that the 
nonlinear matrix elements can be calculated exactly using ap- 
propriately chosen quadrature grids and have described how 
this can be implemented as an efficient algorithm. We have 
also discussed numerous properties of our approach that al- 
low various observables to be calculated, either in the spectral 
basis or by making appropriate use of quadrature grids. Fi- 
nally, we have applied our method to an anisotropic 3D Bose 
gas to indicate typical features of the solutions. 

The rather unique requirement of the PGPE is that it needs 
to be propagated on a small prescribed basis, usually contain- 
ing of order a thousand modes. For this reason the additional 
computational cost of using a spectral method is offset by the 
rather small basis set size. However, recently our method has 
been used to simulate a trapped system with up to 4 x 10 5 
modes 13T1L and is reasonably efficient even for larger prob- 
lems. However, there is much scope for development of the 
algorithm and we hope that the detailed discussion of numer- 
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ics and example application in this paper will stimulate others 
in the community to develop more flexible and general PGPE 
algorithms 

The algorithm outlined in this paper is an important step to 
realizing a comprehensive theory for simulating the dynamics 
of a trapped finite temperature Bose gas. The next steps will 
involve moving beyond the PGPE description to the Stochas- 
tic Gross-Pitaevskii Equation (SGPE) [10], which includes the 
effects of coupling between the classical and incoherent re- 
gions. This theory differs from the PGPE primarily through 
two new features: (i) Stochastic noise terms effect the evolu- 
tion and cause the transfer of particles and energy into and out 
of the classical region. These contributions have already been 
added to the numerical approach outline here [39]. (ii) Scat- 
tering terms corresponding to collisions between incoherent 



and classical region atoms that cause the transfer of momen- 
tum into the system. These terms are more challenging to 
implement and are the subject of current work. 
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